#include "catch.hpp"
#include "inelastic/e_ep_mu_util/ep_mu.h"
#include <range/v3/all.hpp>
#include "generalTools/testing.h"
#include <range/v3/all.hpp>

template <typename RangeInt, typename Range, typename Float>
void checkFirstEprime( const RangeInt& jbetaVec, int lat, const Float& enow, const Range& betas, const Range& x, const Float& tev, const RangeInt& finalJbeta, const Range& finalEp ){
  int jbeta;
  Float ep;
  for ( size_t i = 0; i < jbetaVec.size(); ++i ){
    jbeta = jbetaVec[i];
    ep = findFirstEprime( lat, jbeta, enow, betas, x, tev );
    REQUIRE( finalJbeta[i] == jbeta );
    REQUIRE( finalEp[i] == Approx(ep).epsilon(1e-6) );
  }
}





TEST_CASE( "do we need a midpoint" ){ 
  std::vector<double> x(20,0.0), y(20*65,0.0),
  initialY { 195158.09939, -0.87238772, -0.61817875, -0.365456606, -0.11418564, 
               0.13568062,  0.38420180,  0.63141488,  0.87738226 };
  x[0] = 2.54E-3;
  
  for (size_t i = 0; i < initialY.size(); ++i){ y[i*x.size()+0] = initialY[i]; }



  double xm = 1.27e-3;
  std::vector<double> s { -0.868068168, -0.606773245, -0.349269745, -0.095577908, 
  0.154365216, 0.400543547, 0.6430299311, 0.881812004, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0 };

  int i = 2;
  int nbin = 8;
  double tol = 1.0, 
         pdf = 235586.4814;
  std::vector<double> tolerancesVec { 1.0, 0.01, 0.5, 0.56, 0.59 };
  std::vector<bool>   answersVec { false, true, true, true, false };
  for ( size_t j = 0; j < tolerancesVec.size(); ++j ){
    REQUIRE( needMidpoint(x,y,xm,i,nbin,s,tolerancesVec[j],pdf) == answersVec[j] );
  }

} // TEST CASE



/*
*/


TEST_CASE( "do 330" ){ 


  int imax = 20;
  std::vector<double> x(imax,0.0), y(imax*65,0.0), initialY { 194040.3395, 
                     -0.745374226, -0.240023402, 0.259736457, 0.754304965 };
  x[0] = 2.5607298E-3;
  for (size_t i = 0; i < initialY.size(); ++i){ y[i*x.size()+0] = initialY[i]; }

  int lat = 0, iinc = 2, lasym = 0;
  double jnz = 0.0;
  double tev = 2.5507297688e-2, tol = 5.0E-2;
  double enow;
  std::vector<double> alphas { 1.1, 2.2, 3.3, 4.5, 5.8 },
                       betas { 0.1, 0.2, 1.3, 1.4, 2.5, 2.6, 3.7 },
                       sab { -0.18259619, -0.30201347, -3.93654779, -3.98809174, 
  -4.33545607, -4.39515402, -5.88934921, -0.76225291, -0.81658341, -3.14161459, 
  -3.30566186, -3.90554652, -3.96233362, -5.23696660, -1.19182884, -1.23155471, 
  -2.79610565, -2.95633099, -3.74989225, -3.80837585, -4.93373911, -1.58342860, 
  -1.61310713, -2.71233943, -2.84291608, -3.69699590, -3.75199349, -4.77433858, 
  -1.96121202, -1.98720663, -2.78454600, -2.88531460, -3.71288120, -3.77142141, 
  -4.71158392 };

  double az = 0.99917, sigma_b = 163.72792237, sigma_b2 = 0.0, teff = 0.120441926;
  int nbin = 4, jbeta = 1,j = 0;
  std::vector<double> boundXsVec {sigma_b,sigma_b2};

  std::vector<double> scr(2*imax*65,0.0), xsi(100,0.0);
  


  std::vector<double> correctX(imax), pdf(imax), mu1(imax), mu2(imax), mu3(imax), mu4(imax);
  std::vector<double> lastVals { 0.0, 0.0, 0.0, 0.0, 0.0 };


  GIVEN( "Small initial energy (E=1e-5 eV)" ){
    enow = 1e-5;

    THEN( "Returned x values, y vector, and moment values are correct" ){
      do_330(enow,x,y,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,xsi,lastVals,jnz);

      correctX = { 2.560729E-3, 1.920547E-3, 1.280364E-3, 3.200912E-4, 1.600456E-4, 
      8.002281E-5, 4.001140E-5, 2.000570E-5, 1.000285E-5, 5.001426E-6, 2.500713E-6, 
      1.250356E-6, 6.251782E-7, 3.125891E-7, 1.562945E-7, 7.814728E-8, 3.907364E-8, 
      1.953682E-8, 9.768411E-9, 0.0 };
      pdf = { 194040.3, 213357.4, 234616.8, 270904.3, 278088.3, 281762.2, 
           284435.9, 284849.7, 281681.0, 201641.8, 142771.2, 100595.0, 71050.27, 
           50112.44, 35395.39, 25015.29, 17684.03, 12502.95, 8840.380, 0.0};
      mu1 = {     -0.7453742, -0.7422468, -0.7375231, -0.7158896, -0.6992124, 
      -0.6768329, -0.6449507, -0.6021970, -0.5463158, -0.6019580, -0.6442864, 
      -0.6756062, -0.6972283, -0.7129893, -0.7239285, -0.7315982, -0.7369996, 
      -0.7408188, -0.7435070, 0.00000000 };
      mu2 = {      -0.24002340, -0.23307917, -0.22265027, -0.17496084, -0.13820280, 
      -0.08941394, -0.02005056,  0.07379239,  0.19662745,  0.07431822, -0.01860502, 
      -0.08672420, -0.13381656, -0.16855131, -0.19263496, -0.20951361, -0.22139824, 
      -0.22977689, -0.23570553,  0.0 };
      mu3 = {    0.2597364, 0.2667551, 0.2772268, 0.3250109, 0.3618243, 0.4093513, 
      0.4769398, 0.5740460, 0.6965196, 0.5745724, 0.4783592, 0.4120090, 0.3662356, 
      0.3314537, 0.3073638, 0.2904852, 0.2786010, 0.2702225, 0.2643070, 0.0 };
      mu4 = {    0.7543049, 0.7574966, 0.7623170, 0.7840282, 0.8005880, 0.8224674, 
      0.8530495, 0.8969429, 0.9528127, 0.8971790, 0.8537006, 0.8236852, 0.8025615, 
      0.7869531, 0.7760527, 0.7683948, 0.7629976, 0.7591794, 0.7564790, 0.0 };

      REQUIRE( 20 == j );
      REQUIRE( 2  == jbeta );

      REQUIRE( ranges::equal(x,correctX,equal) );
      std::vector<std::vector<double>> correct { pdf, mu1, mu2, mu3, mu4 };
      for ( size_t i = 0; i < correct.size(); ++i ){
        for ( size_t j = 0; j < correct[i].size(); ++j ){
          REQUIRE( correct[i][j] == Approx(y[i*x.size()+j]).epsilon(1e-6) );
        }
      }
      for ( size_t i = correct.size()*x.size(); i < y.size(); ++i){
        REQUIRE( 0.0 == Approx(y[i]).epsilon(1e-6) );
      }

      std::vector<double> correctSCR { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9.768411E-9, 
      8.8403807E3,-7.435070E-1,-2.357055E-1, 2.643070E-1, 7.564790E-1, 1.953682E-8, 
      1.2502953E4,-7.408188E-1,-2.297769E-1, 2.702225E-1, 7.591794E-1, 3.907364E-8, 
      1.7684034E4,-7.369996E-1,-2.213982E-1, 2.786010E-1, 7.629976E-1, 7.814728E-8, 
      2.5015295E4,-7.315982E-1,-2.095136E-1, 2.904852E-1, 7.683948E-1, 1.562945E-7, 
      3.5395394E4,-7.239285E-1,-1.926349E-1, 3.073638E-1, 7.760527E-1, 3.125891E-7, 
      5.0112448E4,-7.129893E-1,-1.685513E-1, 3.314537E-1, 7.869531E-1, 6.251782E-7, 
      7.1050277E4,-6.972283E-1,-1.338165E-1, 3.662356E-1, 8.025615E-1, 1.250356E-6, 
      1.0059506E5,-6.756062E-1,-8.672421E-2, 4.120090E-1, 8.236852E-1, 2.500713E-6, 
      1.4277126E5,-6.442864E-1,-1.860502E-2, 4.783592E-1, 8.537006E-1, 5.001426E-6, 
      2.0164182E5,-6.019580E-1, 7.431822E-2, 5.745724E-1, 8.971790E-1, 1.000285E-5, 
      2.8168109E5,-5.463158E-1, 1.966274E-1, 6.965196E-1, 9.528127E-1, 2.000570E-5, 
      2.8484972E5,-6.021970E-1, 7.379239E-2, 5.740460E-1, 8.969429E-1, 4.001140E-5, 
      2.8443599E5,-6.449507E-1,-2.005056E-2, 4.769398E-1, 8.530495E-1, 8.002281E-5, 
      2.8176229E5,-6.768329E-1,-8.941394E-2, 4.093513E-1, 8.224674E-1, 1.600456E-4, 
      2.7808831E5,-6.992124E-1,-1.382028E-1, 3.618243E-1, 8.005880E-1, 3.200912E-4, 
      2.7090434E5,-7.158896E-1,-1.749608E-1, 3.250109E-1, 7.840282E-1, 6.401824E-4, 
      2.5808226E5,-7.280949E-1,-2.018525E-1, 2.980850E-1, 7.718158E-1, 1.280364E-3, 
      2.3461686E5,-7.375232E-1,-2.226502E-1, 2.772268E-1, 7.623170E-1, 1.920547E-3, 
      2.1335745E5,-7.422468E-1,-2.330791E-1, 2.667551E-1, 7.574966E-1 };

        for ( size_t i = 0; i < correctSCR.size(); ++i ){
          REQUIRE( scr[i] == Approx(correctSCR[i]).epsilon(1e-6) );
        }

    } // THEN
  } // GIVEN



  GIVEN( "Medium initial energy (E=1e-3 eV)" ){
    enow = 1e-3;

    THEN( "Returned x values, y vector, and moment values are correct" ){

      do_330(enow,x,y,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,xsi,lastVals,jnz);

      correctX = { 2.5607298E-3, 2.5607297E-3, 2.5607295E-3, 6.4018245E-4, 
      1.6004562E-4, 8.0022810E-5, 4.0011405E-5, 2.0005703E-5, 1.0002852E-5, 
      5.0014260E-6, 2.5007130E-6, 1.2503565E-6, 6.2517825E-7, 3.1258913E-7, 
      1.5629457E-7, 7.8147285E-8, 3.9073643E-8, 1.9536822E-8, 9.7684110E-9, 0.0 };
      pdf = { 194040.3, 22928.53, 22928.53, 22004.47, 10421.57, 7271.812, 
      5109.760, 3602.049, 2543.152, 1796.915, 1270.129, 897.9471, 634.8844, 
      448.9098, 317.4197, 224.4469, 158.7070, 112.2224, 79.35316, 0.0 };
      mu1 = {      -0.7453742, -0.6639123, -0.6639123, -0.5966787, -0.6825695, 
      -0.7041852,  -0.7182603, -0.7278000, -0.7343954, -0.7390024, -0.7422431, 
      -0.7445198,  -0.7461270, -0.7472623, -0.7480646, -0.7486316, -0.7490325, 
      -0.7493159,  -0.7495163, 0.0 };
      mu2 = {      -0.24002340, -0.06436183, -0.06436181, 0.08567835, -0.10222948, 
      -0.14962484, -0.18042466, -0.20130111, -0.21574559, -0.22584518, -0.23293962, 
      -0.23794773, -0.24148212, -0.24397868, -0.24574300, -0.24699014, -0.24787184, 
      -0.24849522, -0.24893599, 0.0 };
      mu3 = {    0.2597364, 0.4272364, 0.4272364, 0.5853635, 0.3965669, 0.3496719, 
      0.3192034, 0.2985083, 0.2841580, 0.2741063, 0.2670555, 0.2620485, 0.2585157, 
      0.2560201, 0.2542564, 0.2530096, 0.2521280, 0.2515047, 0.2510640, 0.0 };
      mu4 = {    0.7543049, 0.8233353, 0.8233353, 0.9011784, 0.8153225, 0.7947936, 
      0.7812307, 0.7719456, 0.7654775, 0.7609340, 0.7577044, 0.7554555, 0.7538610, 
      0.7527317, 0.7519323, 0.7513667, 0.7509666, 0.7506835, 0.7504834, 0.0 };

      REQUIRE( 43 == j );
      REQUIRE( 2  == jbeta );

      REQUIRE( ranges::equal(x,correctX,equal) );
      std::vector<std::vector<double>> correct { pdf, mu1, mu2, mu3, mu4 };
      for ( size_t i = 0; i < correct.size(); ++i ){
        for ( size_t j = 0; j < pdf.size(); ++j ){
          REQUIRE( correct[i][j] == Approx(y[i*x.size()+j]).epsilon(1e-6) );
        }
      }
      for ( size_t i = correct.size()*x.size(); i < y.size(); ++i){
        REQUIRE( 0.0 == Approx(y[i]).epsilon(1e-6) );
      }
  
      std::vector<double> correctSCR {0.0, 0.0,  0.0,  0.0,  0.0,  0.0,  9.7684110E-9,  
      7.935317E1, -7.495163E-1, -2.489359E-1,  2.510640E-1,  0.7504834,  1.953682E-8,  
      1.122224E2, -7.493159E-1, -2.484952E-1,  2.515047E-1,  0.7506835,  3.907364E-8,  
      1.587070E2, -7.490325E-1, -2.478718E-1,  2.521281E-1,  0.7509666,  7.814728E-8,  
      2.244469E2, -7.486316E-1, -2.469901E-1,  2.530096E-1,  0.7513667,  1.562945E-7,  
      3.174197E2, -7.480646E-1, -2.457430E-1,  2.542564E-1,  0.7519323,  3.125891E-7,  
      4.489098E2, -7.472623E-1, -2.439786E-1,  2.560201E-1,  0.7527317,  6.251782E-7,  
      6.348844E2, -7.461270E-1, -2.414821E-1,  2.585157E-1,  0.7538610,  1.250356E-6,  
      8.979471E2, -7.445198E-1, -2.379477E-1,  2.620485E-1,  0.7554555,  2.500713E-6,  
      1.270129E3, -7.422431E-1, -2.329396E-1,  2.670555E-1,  0.7577044,  5.001426E-6,  
      1.796915E3, -7.390024E-1, -2.258451E-1,  2.741063E-1,  0.7609340,  1.000285E-5,  
      2.543152E3, -7.343954E-1, -2.157456E-1,  2.841580E-1,  0.7654776,  2.000570E-5,  
      3.602049E3, -7.278000E-1, -2.013011E-1,  2.985083E-1,  0.7719456,  4.001140E-5,  
      5.109760E3, -7.182603E-1, -1.804246E-1,  3.192034E-1,  0.7812307,  8.002281E-5,  
      7.271812E3, -7.041852E-1, -1.496248E-1,  3.496719E-1,  0.7947936,  1.600456E-4,  
      1.042157E4, -6.825695E-1, -1.022294E-1,  3.965669E-1,  0.8153225,  3.200912E-4,  
      1.505175E4, -6.496485E-1, -3.072243E-2,  4.657100E-1,  0.8468844,  6.401824E-4,  
      2.200447E4, -5.966787E-1,  8.567835E-2,  5.853635E-1,  0.9011784,  9.602736E-4,  
      2.786274E4, -5.468695E-1,  1.954414E-1,  6.954181E-1,  0.9523540,  1.280364E-3,  
      2.744101E4, -5.771819E-1,  1.285756E-1,  6.282528E-1,  0.9211990,  1.920547E-3,  
      2.500692E4, -6.313135E-1,  8.429071E-3,  5.050226E-1,  0.8609199,  2.240638E-3,  
      2.378573E4, -6.513582E-1, -3.717092E-2,  0.4529485, 0.8359541,  2.400684E-3,  
      2.346375E4, -6.564133E-1, -4.759721E-2,  0.4444134, 0.8318966,  2.480707E-3,  
      2.319436E4, -6.602429E-1, -5.614899E-2,  0.4356622, 0.8275429,  2.520718E-3,  
      2.306101E4, -6.620969E-1, -6.029621E-2,  0.4314102, 0.8254215,  2.540724E-3,  
      2.299466E4, -6.630093E-1, -6.233904E-2,  0.4293137, 0.8243741,  2.550727E-3,  
      2.296156E4, -6.634620E-1, -6.335292E-2,  0.4282727, 0.8238536,  2.555728E-3,  
      2.294504E4, -6.636874E-1, -6.385800E-2,  0.4277539, 0.8235942,  2.558229E-3,  
      2.293678E4, -6.637999E-1, -6.411007E-2,  0.4274950, 0.8234647,  2.559479E-3,  
      2.293265E4, -6.638561E-1, -6.423600E-2,  0.4273657, 0.8234000,  2.560104E-3,  
      2.293059E4, -6.638842E-1, -6.429893E-2,  0.4273010, 0.8233676,  2.560417E-3,  
      2.292956E4, -6.638983E-1, -6.433039E-2,  0.4272687, 0.8233514,  2.560573E-3,  
      2.292904E4, -6.639053E-1, -6.434612E-2,  0.4272525, 0.8233434,  2.560651E-3,  
      2.292878E4, -6.639088E-1, -6.435398E-2,  0.4272445, 0.8233393,  2.560690E-3,  
      2.292865E4, -6.639105E-1, -6.435791E-2,  0.4272404, 0.8233373,  2.560710E-3,  
      2.292859E4, -6.639114E-1, -6.435988E-2,  0.4272384, 0.8233363,  2.560720E-3,  
      2.292856E4, -6.639119E-1, -6.436086E-2,  0.4272374, 0.8233358,  2.560725E-3,  
      2.292854E4, -6.639121E-1, -6.436136E-2,  0.4272369, 0.8233355,  2.560727E-3,  
      2.292853E4, -6.639122E-1, -6.436160E-2,  0.4272366, 0.8233354,  2.560728E-3,  
      2.292853E4, -6.639122E-1, -6.436172E-2,  0.4272365, 0.8233353,  2.560729E-3,  
      2.292853E4, -6.639123E-1, -6.436178E-2,  0.4272365, 0.8233353,  2.560729E-3,  
      2.292853E4, -6.639123E-1, -6.436181E-2,  0.4272364, 0.8233353,  2.560729E-3,  
      2.292853E4, -6.639123E-1, -6.436183E-2,  0.4272364, 0.8233353,  0.0,  0.0};
        for ( size_t i = 0; i < correctSCR.size(); ++i ){
          REQUIRE( scr[i] == Approx(correctSCR[i]).epsilon(1e-6) );
        }






    } // THEN 
  } // GIVEN

} // TEST CASE


TEST_CASE( "313" ) {
  int lat = 1, jbeta = -7;
  double enow, ep, tev = 2.5507297688E-2;
  std::vector<double> betas { 0.1, 0.2, 1.3, 1.4, 2.5, 2.6, 3.7 }, x(20,0.0);


  GIVEN( "Various jbeta values, fixed enow" ){ 
    std::vector<int> jbetaVec(15), finalJbeta(15);
    std::vector<double> finalEp(15);

    jbetaVec    = ranges::view::iota(-7,8);

    enow = 1e-5;
    finalJbeta = {1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7};
    finalEp    = {2.54e-3, 2.54e-3, 2.54e-3, 2.54e-3, 2.54e-3, 2.54e-3, 2.54e-3, 
                   2.54e-3, 2.54e-3, 5.07e-3, 3.29e-2, 3.543e-2, 6.326e-2, 
                   6.579e-2, 9.362e-2};
    checkFirstEprime( jbetaVec, lat, enow, betas, x, tev, finalJbeta, finalEp );


    enow = 1e-3;
    finalJbeta = {1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7};
    finalEp    = {3.53e-3, 3.53e-3, 3.53e-3, 3.53e-3, 3.53e-3, 3.53e-3, 3.53e-3, 
                  3.53e-3, 3.53e-3, 6.06e-3, 3.389e-2, 3.642e-2, 6.425e-2, 
                  6.678e-2, 9.461e-2};
    checkFirstEprime( jbetaVec, lat, enow, betas, x, tev, finalJbeta, finalEp );


    enow = 1e-1;
    finalJbeta = {-7, -6, -5, -4, -3, -2, -1, 1, 1, 2, 3, 4, 5, 6, 7 };
    finalEp    = {6.39000153e-3, 3.4220001e-2, 3.6750001e-2, 6.4580001e-2, 
                  6.7110001e-2, 9.4940001e-2, 9.7470001e-2, 0.10253,  0.10253, 
                  0.10506, 0.13289, 0.13542, 0.16325, 0.16578, 0.19361 };
    checkFirstEprime( jbetaVec, lat, enow, betas, x, tev, finalJbeta, finalEp );

  } // GIVEN



} // TEST CASE







TEST_CASE( "do 330 (and some things around it)" ){ 
  int imax = 20, lat = 0, iinc = 2, lasym = 0;
  double tev = 2.5507297688e-2, tol = 5.0E-2;
  std::vector<double> 
  alphas { 1.1, 2.2, 3.3, 4.5, 5.8 },
  betas { 0.1, 0.2, 1.3, 1.4, 2.5, 2.6, 3.7 },
  sab {-0.18259619, -0.30201347, -3.93654779, -3.98809174, -4.33545607, 
  -4.39515402, -5.88934921, -0.76225291, -0.81658341, -3.14161459, -3.30566186, 
  -3.90554652, -3.96233362, -5.23696660, -1.19182884, -1.23155471, -2.79610565, 
  -2.95633099, -3.74989225, -3.80837585, -4.93373911, -1.58342860, -1.61310713, 
  -2.71233943, -2.84291608, -3.69699590, -3.75199349, -4.77433858, -1.96121202, 
  -1.98720663, -2.78454600, -2.88531460, -3.71288120, -3.77142141, -4.71158392 };

  std::vector<double> y(20*65,0.0);


  double az = 0.99917, sigma_b = 163.72792237, sigma_b2 = 0.0, teff = 0.120441926;
  int nbin = 4, jbeta = 1,j = 0;
  std::vector<double> boundXsVec {sigma_b,sigma_b2};


  std::vector<double> scr(65*imax*8,0.0);
  std::vector<double> xsi(4,0.0), correctVals(4);

  std::vector<double>  correctSCR;
  std::vector<double> lastVals { 0.0, 0.0, 0.0, 0.0, 0.0 };


  GIVEN( "2 bins requested" ){
    double enow,tol;
    nbin = 2;



    WHEN( "initial energy is small (E = 1e-6 eV)" ){ 
      enow = 1e-6; 
      AND_WHEN( "various tolerances considered" ){ 
        tol = 5e-1;
        for ( auto& val : scr ){ val = 0.0; }
        auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
        correctSCR = { 0, 0, 0, 0, 1.245962E-6, 134.9500, -0.187640, 0.806634, 
        2.491924E-6, 134.7343,-0.275233,  0.72264533, 4.98384E-6, 133.6925, 
       -0.3407188,   0.658092, 9.967695E-6, 131.7394,-0.391743,   0.608060, 
        1.993539E-5, 130.9234,-0.424779,  0.57516619, 3.98707E-5, 130.2870, 
       -0.4474311,   0.552549, 7.974156E-5, 129.4080,-0.463325,   0.536666, 
        1.594831E-4, 127.8354,-0.474647,  0.52534821, 3.18966E-4, 124.8242, 
       -0.4828510,   0.517144, 6.379325E-4, 119.0482,-0.488959,   0.511033, 
        1.275865E-3, 108.3016,-0.493721,  0.50626439, 2.55173E-3, 89.63757, 
       -0.4977176,   0.502254, 3.827095E-3, 74.19671,-0.499894,   0.500063, 
        5.102460E-3, 61.41592,-0.501429,  0.49851398, 1.21169E-2, 19.07564, 
       -0.5017302,   0.498257, 1.913147E-2, 8.602442,-0.502753,   0.497225, 
        3.316049E-2, 2.080990,-0.504043,  0.49592504, 3.57112E-2, 2.032783, 
       -0.5033939,   0.496583, 6.376924E-2, 1.440228,-0.501071,   0.498924, 
        6.631997E-2, 1.339506,-0.501021,  0.49897454, 9.43780E-2, 0, 0, 0 };
        for ( size_t i = 0; i < correctSCR.size(); ++i ){
          REQUIRE( scr[i] == Approx(correctSCR[i]).epsilon(1e-6) );
        }
        correctVals = {6843.981, 2.3013E-3, -1.2485E-1,  8.9429E-4 };
        REQUIRE( ranges::equal( correctVals, out, equal ) );

        tol = 1.0; 
        jbeta = -7; j = 0;
        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : xsi ){ val = 0.0; }
        for ( auto& val : lastVals){ val = 0.0; }
        out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
        correctSCR = {0, 0, 0, 0, 2.551730E-3, 94.58360, -0.4977176, 0.5022543, 
        3.827095E-3, 78.29074, -0.4998946, 0.50006329, 5.102460E-3, 64.80473, 
        -0.5014298, 0.49851398, 0.01913147, 9.077108, -0.5027530, 0.49722554, 
        0.03316049, 2.195815, -0.5040439, 0.49592504, 0.03571122, 2.144948, 
        -0.5033939, 0.49658361, 0.06376924, 1.519697, -0.5010712, 0.49892435, 
        0.06631997, 1.413417, -0.5010211, 0.49897454, 0.09437800, 0, 0, 0 };
        correctVals = { 6486.092, -8.8873e-4, -0.12503, -3.3313e-4 };
        REQUIRE( ranges::equal( correctSCR,  scr, equal ) );
        REQUIRE( ranges::equal( correctVals, out, equal ) );

      } // AND WHEN
    } // WHEN

    WHEN( "initial energy is high (E = 5eV)" ){
      enow = 5.0;
      AND_WHEN( "various tolerances considered" ){
        tol = 2.0;
        jbeta = -7;
        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : xsi ){ val = 0.0; }
        for ( auto& val : lastVals ){ val = 0.0; }
        auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
        correctSCR = {0.0, 0.0, 0.0, 0.0, 4.90562, 0.373818, 0.925628, 0.983460, 
        4.933680, 0.433181, 0.935837, 0.987965, 4.936231, 0.425528, 0.935521, 
        0.987938, 4.964289, 0.288884, 0.925321, 0.982937, 4.966840, 0.232108, 
        0.917665, 0.982192, 4.994898, 1.459286, 0.975571, 0.997333, 4.997449, 
        1.787917, 0.979322, 0.998269, 5.002550, 1.617792, 0.979333, 0.998270,
        5.005101, 1.194692, 0.975595, 0.997336, 5.033159, 0.063227, 0.918258, 
        0.982317, 5.035710, 0.071134, 0.925829, 0.983050, 5.063768, 0.034977, 
        0.936428, 0.988081, 5.066319, 0.032224, 0.936778, 0.988112,5.094377,0,0,0};
        correctVals = { 16.76715, 0.95637, 0.87324, 0.75837 };
        REQUIRE( ranges::equal(correctSCR,  scr, equal) );
        REQUIRE( ranges::equal(correctVals, out, equal) );


        j = 0;
        tol = 1e-1;
        jbeta = -7;
        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : xsi ){ val = 0.0; }
        out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
        correctSCR = { 0.0, 0.0, 0.0, 0.0, 1.87134E-5, 2.94462E-3, -0.499319, 
        0.5006542, 3.74269E-5, 4.16408E-3, -0.499030, 0.500917, 7.48538E-5, 
        5.88825E-3, -0.498612, 0.501284, 1.49707E-4, 8.32540E-3, -0.498010, 
        0.501781, 2.99415E-4, 0.0117688, -0.497120, 0.502466, 5.98831E-4, 
        0.0166291, -0.495808, 0.503360, 1.19766E-3, 0.0234759, -0.493836, 
        0.5044974, 2.395324E-3, 0.0330834, -0.4908072, 0.505861, 4.790648E-3, 
        0.0464594, -0.4860256, 0.507323, 9.58129E-3, 0.0647928, -0.478226, 
        0.508522, 0.0191625, 0.089150, -0.465005, 0.508693, 0.0383251, 0.119584, 
        -0.441645, 0.506586, 0.076650, 0.153343, -0.398784, 0.501180, 0.1533007, 
        0.184997, -0.315236, 0.496817, 0.3066014, 0.202430, -0.171198, 0.497999, 
        0.613202, 0.207263, 0.036639, 0.523631, 1.226405, 0.208707, 0.283923, 
        0.606369, 2.452811, 0.207073, 0.569175, 0.760800, 4.905623, 0.174472, 
        0.933450, 0.985000, 4.933681, 0.199517, 0.940955, 0.989624, 4.936231, 
        0.178730, 0.936943, 0.988202, 4.964289, 0.202265, 0.947704, 0.990366, 
        4.966840, 0.209896, 0.949741, 0.990666, 4.980869, 0.306808, 0.961109, 
        0.993929, 4.987884, 0.416525, 0.969360, 0.995829, 4.994898, 0.646950, 
        0.976223, 0.997509, 4.997449, 0.753987, 0.980776, 0.998317, 4.998724, 
        0.793903, 0.981686, 0.998638, 4.999362, 0.813694, 0.982116, 0.998789, 
        4.999681, 0.819510, 0.982256, 0.998841, 4.999840, 0.820134, 0.982285, 
        0.998856, 4.999920, 0.819731, 0.982286, 0.998859, 4.999960, 0.819458, 
        0.982286, 0.998861, 4.999970, 0.834126, 0.980848, 0.998851, 4.999975, 
        0.853587, 0.981134, 0.998887, 4.999977, 0.866309, 0.981363, 0.998945, 
        4.999978, 0.874852, 0.981514, 0.998983, 4.999979, 0.944516, 0.980657, 
        0.998841, 4.999979, 1.082571, 0.977254, 0.998467, 4.999979, 1.370069, 
        0.966263, 0.997361, 4.999980, 1.370998, 0.966280, 0.997369, 4.999980, 
        1.976504, 0.937578, 0.993046, 4.999980, 1.977461, 0.937596, 0.993061, 
        4.999980, 5.743719, 0.742946, 0.947992, 4.999980, 5.745691, 0.742990, 
        0.948056, 4.999980, 5.748759, 0.743057, 0.948155, 4.999981, 5.755301, 
        0.743201, 0.948365, 4.999982, 5.770053, 0.743525, 0.948835, 4.999985, 
        5.808005, 0.744356, 0.950028, 4.999990, 5.906783, 0.746501, 0.953016, 
        5.000000, 6.042281, 0.749423, 0.956880, 5.000010, 5.902344, 0.746455, 
        0.952953, 5.000015, 5.802889, 0.744319, 0.949975, 5.000017, 5.764778, 
        0.743497, 0.948794, 5.000018, 5.748816, 0.743152, 0.948293, 5.000019, 
        5.742325, 0.743012, 0.948088, 5.000019, 5.739280, 0.742946, 0.947992, 
        5.000019, 1.975926, 0.937596, 0.993061, 5.000019, 1.974962, 0.937578, 
        0.993046, 5.000020, 1.369924, 0.966280, 0.997369, 5.000020, 0.942044, 
        0.980629, 0.998833, 5.000021, 0.872575, 0.981487, 0.998976, 5.000022, 
        0.864291, 0.981341, 0.998940, 5.000025, 0.851985, 0.981120, 0.998883, 
        5.000030, 0.826746, 0.980726, 0.998816, 5.000039, 0.818175, 0.982286, 
        0.998861, 5.000079, 0.817173, 0.982287, 0.998859, 5.000159, 0.815026, 
        0.982285, 0.998856, 5.000318, 0.809335, 0.982257, 0.998841, 5.000637, 
        0.793556, 0.982119, 0.998789, 5.001275, 0.755153, 0.981691, 0.998638, 
        5.002550, 0.682243, 0.980787, 0.998318, 5.005101, 0.529706, 0.976249, 
        0.997511, 5.012116, 0.259032, 0.969439, 0.995840, 5.019130, 0.144898, 
        0.961263, 0.993954, 5.026145, 0.088623, 0.954869, 0.992231, 5.033159, 
        0.057206, 0.950075, 0.990727, 5.035710, 0.049868, 0.948076, 0.990433, 
        5.049739, 0.028259, 0.943855, 0.990186, 5.056753, 0.018185, 0.936482, 
        0.988095, 5.063768, 0.016338, 0.942060, 0.989837, 5.066319, 0.014814, 
        0.941800, 0.989760, 5.080348, 8.194329E-3, 0.938315, 0.988339, 5.087362, 
        5.518491E-3, 0.9326195, 0.98594691, 5.094377, 0.0, 0.0, 0.0};
        correctVals = { 39.08391, 0.61572, 0.20856, 0.016571 };
        REQUIRE( ranges::equal(correctSCR,  scr, equal) );
        REQUIRE( ranges::equal( correctVals, out, equal ) );
      } // AND WHEN
    } // WHEN




    WHEN( "higher energy, high tolerance" ){
      enow = 0.5;
      tol = 5e-1;
      jbeta = -7;
      for ( auto& val : scr ){ val = 0.0; }

      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

      correctSCR = {0.0, 0.0, 0.0, 0.0, 0.2028115, 2.088333, -0.0178393, 
      0.720697, 0.405623, 1.684778, 0.273684, 0.827656, 0.433681, 1.922670, 
      0.374357, 0.877023, 0.436231, 1.894687, 0.373533, 0.876688, 0.464289, 
      2.185131, 0.474614, 0.900139, 0.466840, 2.254948, 0.490033, 0.902993, 
      0.494898, 6.813146, 0.779380, 0.974165, 0.497449, 8.425869, 0.814161, 
      0.982932, 0.498724, 8.878005, 0.823191, 0.986081, 0.499362, 9.109986, 
      0.827574, 0.987619, 0.499681, 9.227942, 0.829741, 0.988389, 0.499840, 
      9.277499, 0.830671, 0.988728, 0.499920, 9.288740, 0.830934, 0.988836, 
      0.499960, 9.288837, 0.830984, 0.988865, 0.500000, 34.13863, 0.721390, 
      0.938708, 0.500039, 9.274387, 0.830995, 0.988865, 0.500079, 9.259794, 
      0.830956, 0.988837, 0.500159, 9.219706, 0.830715, 0.988731, 0.500318, 
      9.113413, 0.829829, 0.988395, 0.500637, 8.885335, 0.827754, 0.987633, 
      0.501275, 8.445616, 0.823558, 0.986113, 0.502550, 7.624984, 0.814930, 
      0.983011, 0.505101, 5.583066, 0.781296, 0.974422, 0.519130, 1.526250, 
      0.632802, 0.939442, 0.533159, 0.612663, 0.521724, 0.909378, 0.535710, 
      0.538822, 0.503196, 0.906930, 0.563768, 0.157392, 0.436427, 0.890570, 
      0.566318, 0.144765, 0.439953, 0.891308, 0.594377, 0.0, 0.0, 0.0};

      correctVals = { 36.29414, 0.56012, 0.16052, 0.060581 };
      REQUIRE( ranges::equal( correctVals, out, equal ) );
      REQUIRE( ranges::equal( correctSCR,  scr, equal ) );

    } // WHEN

    WHEN( "rather high energy (E = 2 eV) and rather high tolerance (0.5)" ){
      enow = 2.0;
      tol = 5e-1;
      jbeta = -7;
      for ( auto& val : scr ){ val = 0.0; }
      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
      correctSCR = {0.0, 0.0, 0.0, 0.0, 0.952811, 0.660070, 0.415912, 0.761437, 
      1.905623, 0.520008, 0.809530, 0.958128, 1.933681, 0.585289, 0.833890, 
      0.969934, 1.936231, 0.576276, 0.833329, 0.969824, 1.964289, 0.658464, 
      0.858282, 0.975389, 1.966840, 0.679192, 0.862343, 0.976092, 1.994898, 
      2.087298, 0.939601, 0.993287, 1.997449, 2.510947, 0.948097, 0.995539, 
      1.998724, 2.634583, 0.950412, 0.996384, 1.999362, 2.701458, 0.951593, 
      0.996801, 1.999681, 2.730235, 0.952102, 0.996982, 1.999840, 2.736710, 
      0.952243, 0.997039, 1.999920, 2.736691, 0.952266, 0.997054, 1.999960, 
      2.735888, 0.952266, 0.997058, 2.000000, 18.67175, 0.739671, 0.943746, 
      2.000039, 2.731618, 0.952267, 0.997058, 2.000079, 2.728161, 0.952268, 
      0.997054, 2.000159, 2.719684, 0.952247, 0.997039, 2.000318, 2.696297, 
      0.952110, 0.996982, 2.000637, 2.634811, 0.951608, 0.996801, 2.001275, 
      2.506125, 0.950443, 0.996386, 2.002550, 2.272072, 0.948163, 0.995545, 
      2.005101, 1.708756, 0.939747, 0.993302, 2.019130, 0.457752, 0.898796, 
      0.984534, 2.033159, 0.184982, 0.864700, 0.976495, 2.035710, 0.162256, 
      0.860892, 0.975833, 2.063768, 0.0473657,0.839053, 0.970762, 2.066319, 
      0.0435444, 0.839847, 0.970900, 2.094377, 0.0, 0.0, 0.0 };
      correctVals = { 30.44613, 0.70524, 0.31180, 0.042036 };
      REQUIRE( ranges::equal( correctVals, out, equal ) );
      REQUIRE( ranges::equal( correctSCR,  scr, equal ) );
    } // WHEN
  } // GIVEN 




  GIVEN( "four bins requested" ){
    WHEN( "Small energy" ){
      double enow = 1e-6;

      AND_WHEN( "tolerance is small" ){

        double tol = 5e-1;
        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : lastVals ){ val = 0.0; }
        auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

        correctSCR = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.24596E-6, 134.9500, -0.553199,  
        0.177917, 0.672898, 0.940370, 2.49192E-6, 134.7343, -0.609136,  5.866965E-2, 
        0.557295, 0.887995, 4.98384E-6, 133.6925, -0.650658, -3.077921E-2, 0.469878, 
        0.846307, 9.96769E-6, 131.7394, -0.682402, -0.101083, 0.399095, 0.817025, 
        1.99353E-5, 130.9234, -0.703005, -0.146554, 0.353472, 0.796859, 3.98707E-5, 
        130.2870, -0.717149, -0.177712, 0.322288, 0.782811, 7.97415E-5, 129.4080, 
        -0.727080, -0.199571, 0.300426, 0.772906, 1.59483E-4, 127.8354, -0.734154, 
        -0.215139, 0.284857, 0.765838, 3.18966E-4, 124.8242, -0.739281, -0.226420, 
        0.273575, 0.760713, 6.37932E-4, 119.0482, -0.743101, -0.234816, 0.265191, 
        0.756874, 1.27586E-3, 108.3016, -0.746073, -0.241369, 0.258622, 0.753906, 
        2.55173E-3, 89.63757, -0.748566, -0.246868, 0.253107, 0.751401, 3.82709E-3, 
        74.19671, -0.749923, -0.249865, 0.250097, 0.750028, 5.10246E-3, 61.41592, 
        -0.750880, -0.251979, 0.247971, 0.749056, 1.21169E-2, 19.07564, -0.751078, 
        -0.252382, 0.247606, 0.748907, 1.91314E-2, 8.602442, -0.751715, -0.253790, 
        0.246190, 0.748260, 3.31604E-2, 2.080990, -0.752519, -0.255568, 0.244404, 
        0.747445, 3.57112E-2, 2.032783, -0.752115, -0.254672, 0.245308, 0.747859, 
        6.37692E-2, 1.440228, -0.750668, -0.251474, 0.248522, 0.749326, 6.63199E-2, 
        1.339506, -0.750637, -0.251405, 0.248590, 0.749358, 9.43780E-2,0,0,0,0,0};
        correctVals = { 6843.981, 2.3013E-3, -3.1075E-2, 3.6470E-4 };
        REQUIRE( ranges::equal( correctVals, out, equal ) );
        REQUIRE( ranges::equal( correctSCR,  scr, equal ) );


        tol = 1.5;
        jbeta = -7; j = 0;

        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : lastVals ){ val = 0.0; }
        for ( auto& val : xsi ){ val = 0.0; }
        out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

        correctSCR = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.551730E-3, 94.58360, 
       -0.7485668, -0.2468684, 0.25310755, 0.75140105, 3.827095E-3, 78.29074, 
       -0.7499239, -0.2498653, 0.25009783, 0.75002876, 5.102460E-3, 64.80473, 
       -0.7508800, -0.2519797, 0.24797102, 0.74905694, 0.01913147, 9.077108, 
       -0.7517153, -0.2537907, 0.24619046, 0.74826062, 0.03316049, 2.195815, 
       -0.7525198, -0.2555681, 0.24440473, 0.74744536, 0.03571122, 2.144948, 
       -0.7521157, -0.2546722, 0.24530807, 0.74785916, 0.06376924, 1.519697, 
       -0.7506683, -0.2514740, 0.24852203, 0.74932668, 0.06631997, 1.413417, 
       -0.7506371, -0.2514052, 0.24859099, 0.74935809, 0.09437800, 0, 0, 0, 0, 0 };
        correctVals = { 6486.092, -8.8873E-4, -3.1277E-2, -1.2541E-4 };
        REQUIRE( ranges::equal( correctSCR,  scr, equal ) );
        REQUIRE( ranges::equal( correctVals, out, equal ) );

        tol = 5.0;
        jbeta = -7; j = 0;

        for ( auto& val : scr ){ val = 0.0; }
        for ( auto& val : lastVals ){ val = 0.0; }
        for ( auto& val : xsi ){ val = 0.0; }
        out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

        correctSCR = { 0, 0, 0, 0, 0, 0, 2.551730E-3, 70.35257, -0.7485668, 
        -0.246868, 0.2531075, 0.7514010, 5.10246E-3, 48.20264, -0.7508800, 
        -0.251979, 0.2479710, 0.7490569, 3.31604E-2, 1.633277, -0.7525198, 
        -0.255568, 0.2444047, 0.7474453, 3.57112E-2, 1.595441, -0.7521157, 
        -0.254672, 0.2453080, 0.7478591, 6.37692E-2, 1.130372, -0.7506683, 
        -0.251474, 0.2485220, 0.7493266, 6.63199E-2, 1.051319, -0.7506371, 
        -0.251405, 0.2485909, 0.7493580, 9.43780E-2, 0, 0, 0, 0, 0 };
        correctVals = { 8720.050, -8.8804E-4, -3.1283E-2, -1.2535E-4 };
        REQUIRE( ranges::equal( correctVals, out, equal ) );
        REQUIRE( ranges::equal( correctSCR,  scr, equal ) );



        
      } // AND WHEN
    } // WHEN

    WHEN( "Initial energy E = 1e-5 eV" ){
      double enow = 1e-5;
      tol = 5e-2;
      jbeta = -7;
      j = 0;
      for ( auto& val : scr ){ val = 0.0; }
      for ( auto& val : xsi ){ val = 0.0; }
      for ( auto& val : lastVals ){ val = 0.0; }
      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);
      correctSCR = { 0, 0, 0, 0, 0, 0, 9.768411E-9, 
      4.428672, -0.7435070, -0.2357055, 0.2643070, 0.7564790, 1.95368E-8, 
      6.263472, -0.7408188, -0.2297769, 0.2702225, 0.7591794, 3.90736E-8, 
      8.858983, -0.7369996, -0.2213982, 0.2786010, 0.7629976, 7.81472E-8, 
      12.53165, -0.7315983, -0.2095136, 0.2904852, 0.7683948, 1.56294E-7, 
      17.73166, -0.7239286, -0.1926350, 0.3073638, 0.7760527, 3.12589E-7, 
      25.10430, -0.7129893, -0.1685513, 0.3314537, 0.7869531, 6.25178E-7, 
      35.59330, -0.6972284, -0.1338166, 0.3662356, 0.8025615, 1.25035E-6, 
      50.39404, -0.6756062, -0.0867242, 0.4120090, 0.8236852, 2.50071E-6, 
      71.52261, -0.6442864, -0.0186050, 0.4783592, 0.8537006, 5.00142E-6, 
      101.0144, -0.6019581,  0.0743182, 0.5745724, 0.8971790, 1.00028E-5, 
      141.1108, -0.5463158,  0.1966274, 0.6965196, 0.9528127, 2.00057E-5, 
      142.6981, -0.6021971,  0.0737923, 0.5740460, 0.8969429, 4.00114E-5, 
      142.4909, -0.6449508, -0.0200505, 0.4769398, 0.8530495, 8.00228E-5, 
      141.1515, -0.6768329, -0.0894139, 0.4093513, 0.8224674, 1.60045E-4, 
      139.3110, -0.6992125, -0.1382028, 0.3618243, 0.8005880, 3.20091E-4, 
      135.7121, -0.7158896, -0.1749608, 0.3250109, 0.7840282, 6.40182E-4, 
      129.2887, -0.7280950, -0.2018526, 0.2980850, 0.7718158, 1.28036E-3, 
      117.5335, -0.7375232, -0.2226503, 0.2772268, 0.7623170, 1.92054E-3, 
      106.8834, -0.7422468, -0.2330792, 0.2667551, 0.7574966, 2.56073E-3, 
      97.20633, -0.7453742, -0.2400234, 0.2597364, 0.7543049, 3.19841E-3, 
      88.44101, -0.7477164, -0.2452313, 0.2544634, 0.7518883, 3.83609E-3, 
      80.46737, -0.7496037, -0.2494359, 0.2501968, 0.7499233, 5.11146E-3, 
      66.61351, -0.7525702, -0.2560654, 0.2434452, 0.7467990, 5.11156E-3, 
      48.63349, -0.7500001, -0.2500002, 0.2499998, 0.7499999, 5.11167E-3, 
      48.63286, -0.7500002, -0.2500004, 0.2499996, 0.7499998, 5.11188E-3, 
      48.63160, -0.7500003, -0.2500007, 0.2499992, 0.7499996, 5.11231E-3, 
      48.62909, -0.7500006, -0.2500014, 0.2499985, 0.7499993, 5.11317E-3, 
      48.62405, -0.7500013, -0.2500028, 0.2499971, 0.7499986, 5.11488E-3, 
      48.61399, -0.7500025, -0.2500056, 0.2499942, 0.7499973, 5.11831E-3, 
      48.59387, -0.7500051, -0.2500113, 0.2499885, 0.7499947, 5.12516E-3, 
      48.55364, -0.7500101, -0.2500225, 0.2499771, 0.7499894, 5.13886E-3, 
      48.47326, -0.7500202, -0.2500450, 0.2499544, 0.7499790, 5.16626E-3, 
      48.31283, -0.7500403, -0.2500897, 0.2499090, 0.7499581, 5.22106E-3, 
      47.99323, -0.7500801, -0.2501784, 0.2498190, 0.7499166, 5.33066E-3, 
      47.35917, -0.7501585, -0.2503531, 0.2496418, 0.7498351, 5.54986E-3, 
      46.11191, -0.7503107, -0.2506922, 0.2492980, 0.7496768, 5.98827E-3, 
      43.70278, -0.7505983, -0.2513328, 0.2486484, 0.7493777, 6.86508E-3, 
      39.23052, -0.7511177, -0.2524895, 0.2474755, 0.7488377, 8.61871E-3, 
      31.60939, -0.7519948, -0.2544435, 0.2454940, 0.7479251, 0.01212597, 
      20.69604, -0.7533613, -0.2574904, 0.2424012, 0.7464993, 0.01563322, 
      13.77782, -0.7544366, -0.2598920, 0.2399589, 0.7453716, 0.01914047, 
      9.333232, -0.7553412, -0.2619156, 0.2378971, 0.7444178, 0.02264773, 
      6.427451, -0.7561320, -0.2636877, 0.2360879, 0.7435793, 0.02615498, 
      4.494250, -0.7568408, -0.2652788, 0.2344606, 0.7428238, 0.02966223, 
      3.156512, -0.7563538, -0.2641408, 0.2356766, 0.7434113, 0.03316949, 
      2.258159, -0.7578482, -0.2675061, 0.2322225, 0.7418026, 0.03572022, 
      2.205466, -0.7566034, -0.2647018, 0.2351020, 0.7431442, 0.04974923, 
      2.144491, -0.7542702, -0.2594970, 0.2403881, 0.7455822, 0.05676374, 
      1.880905, -0.7526384, -0.2559083, 0.2439006, 0.7470732, 0.06377824, 
      1.561732, -0.7520956, -0.2546458, 0.2453147, 0.7478539, 0.06632897, 
      1.452495, -0.7519976, -0.2544288, 0.2455333, 0.7479539, 0.08035799, 
      0.825841, -0.7520148, -0.2544723, 0.2454837, 0.7479288, 0.08737250, 
      0.572034, -0.7520276, -0.2544978, 0.2454610, 0.7479197, 0.09438700, 
      0, 0, 0, 0, 0 };

      correctVals = { 1996.170, 8.2011E-3, -2.9984E-2, 1.5592E-3 };
      REQUIRE( ranges::equal( correctVals, out, equal ) );
      REQUIRE( ranges::equal( correctSCR,  scr, equal ) );
    } // WHEN



    WHEN( "Initial energy E = 1e-3 eV" ){
      double enow = 1e-3;
      for ( auto& val : scr ){ val = 0.0; }
      tol = 5e-2;

      jbeta = 1;
      j = 0;
      std::vector<double> scr(100,0.0);
      std::vector<double> xsi(100,0.0);
      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

      correctSCR ={ 0, 0, 0, 0, 0, 0, 1.35449E-8, 0.413853,-0.749430,-0.2487471, 
      0.251252, 0.750569, 2.70899E-8, 0.585278,-0.749194,-0.2482280, 0.251771, 
      0.750804, 5.41798E-8, 0.827712,-0.748860,-0.2474939, 0.252505, 0.751138, 
      1.08359E-7, 1.170572,-0.748388,-0.2464556, 0.253544, 0.751609, 2.16719E-7, 
      1.655465,-0.747720,-0.2449869, 0.255012, 0.752275, 4.33438E-7, 2.341258, 
     -0.746775,-0.2429088, 0.257089, 0.753215, 8.66877E-7, 3.311257,-0.745438, 
     -0.2399678, 0.260029, 0.754544, 1.73375E-6, 4.683439,-0.743544,-0.2358027, 
      0.264192, 0.756420, 3.46751E-6, 6.625122,-0.740860,-0.2298976, 0.270068, 
      0.759094, 6.93502E-6, 9.374266,-0.737032,-0.2215258, 0.278407, 0.762879, 
      1.38700E-5, 13.27120,-0.731581,-0.2095807, 0.290286, 0.768242, 2.77400E-5, 
      18.80817,-0.723745,-0.1924268, 0.307311, 0.775901, 5.54801E-5, 26.71330, 
     -0.712323,-0.1674356, 0.332059, 0.786970, 1.10960E-4, 38.11516,-0.695213, 
     -0.1299720, 0.369102, 0.803361, 2.21920E-4, 54.69209,-0.669624,-0.0743532, 
      0.422876, 0.827609, 4.43841E-4, 79.43539,-0.628670,  0.0153657, 0.514483, 
      0.868525, 8.87682E-4, 117.7772,-0.558048,  0.1707997, 0.670740, 0.941121, 
      1.77536E-3, 113.2958,-0.620794,  0.0319876, 0.529339, 0.872838, 2.66304E-3, 
      100.0673,-0.668388,-0.0744072, 0.416899, 0.818150, 3.55073E-3, 88.16909, 
     -0.699620,-0.1453324, 0.342632, 0.780048, 4.82609E-3, 73.63417,-0.730208, 
     -0.2164941, 0.264780, 0.737776, 5.14493E-3, 70.42146,-0.736215,-0.2306806, 
      0.248839, 0.728928, 5.46377E-3, 67.35297,-0.741769,-0.2438752, 0.233842, 
      0.720522, 6.10146E-3, 61.62057,-0.751746,-0.2677746, 0.206299, 0.704774, 
      6.10156E-3, 46.97668,-0.750000,-0.2500015, 0.249998, 0.749999, 6.10167E-3, 
      46.97604,-0.750001,-0.2500030, 0.249996, 0.749998, 6.10188E-3, 46.97475, 
     -0.750002,-0.2500060, 0.249993, 0.749996, 6.10231E-3, 46.97217,-0.750005, 
     -0.2500120, 0.249986, 0.749992, 6.10317E-3, 46.96703,-0.750009,-0.2500239, 
      0.249972, 0.749985, 6.10488E-3, 46.95673,-0.750019,-0.2500479, 0.249945, 
      0.749971, 6.10831E-3, 46.93615,-0.750039,-0.2500957, 0.249891, 0.749942, 
      6.11516E-3, 46.89502,-0.750079,-0.2501914, 0.249782, 0.749885, 6.12886E-3, 
      46.81286,-0.750158,-0.2503823, 0.249565, 0.749771, 6.15626E-3, 46.64898, 
     -0.750315,-0.2507628, 0.249133, 0.749544, 6.21106E-3, 46.32294,-0.750627, 
     -0.2515184, 0.248275, 0.749093, 6.32066E-3, 45.67772,-0.751244,-0.2530089, 
      0.246583, 0.748203, 6.53986E-3, 44.41444,-0.752444,-0.2559105, 0.243290, 
      0.746471, 6.97827E-3, 41.99430,-0.754723,-0.2614241, 0.237035, 0.743183, 
      7.85508E-3, 37.55877,-0.758875,-0.2714704, 0.225636, 0.737187, 9.60871E-3, 
      30.12541,-0.765948,-0.2886299, 0.206120, 0.726890, 0.0131159, 19.65693, 
     -0.7769744,-0.3155805, 0.175180, 0.710385, 0.0166232, 13.08528,-0.7855291, 
     -0.3367411, 0.150489, 0.696953, 0.0201304, 8.844176,-0.7895341,-0.3501393, 
      0.132399, 0.686408, 0.0236377, 6.048807,-0.7918065,-0.3572893, 0.120984, 
      0.678726, 0.0271449, 4.207837,-0.7966820,-0.3661695, 0.110176, 0.671601, 
      0.0306522, 2.988024,-0.8045808,-0.3843479, 0.0931337,0.663254, 0.0341594, 
      2.181687,-0.8132304,-0.4063657, 0.0676373,0.650060, 0.0367102, 2.095439, 
     -0.8046342,-0.3832756, 0.0979029,0.669578, 0.0507392, 1.951280,-0.7829508, 
     -0.3309273, 0.157654, 0.702195, 0.0647682, 1.398690,-0.7686586,-0.2947295, 
      0.199990, 0.723730, 0.0673189, 1.300457,-0.7674876,-0.2918944, 0.203477, 
      0.726032, 0.0813479, 0.737226,-0.7664191,-0.2895516, 0.205975, 0.727899, 
      0.0883625, 0.509961,-0.7678017,-0.292903, 0.202052, 0.725564, 0.095377,0,0,0,0,0};
      correctVals = { 225.7848, 4.6022E-2, -1.3691E-2, 1.5416E-2 };
      REQUIRE( ranges::equal( correctVals, out, equal ) );
      REQUIRE( ranges::equal( correctSCR, scr, equal ) );
    } // WHEN


    WHEN( "Higher initial energy (0.1 eV)" ){
      double enow = 1e-1;
      tol = 1e-1;
      jbeta = -7;
      for ( auto& val : scr ){ val = 0.0; }

      std::vector<double> xsi(100,0.0);
      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

      correctSCR = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.864013E-7, 2.036626E-2, 
      -7.494533E-1, -0.2487979, 0.251201, 0.750545, 1.37280E-6, 2.88021E-2, 
      -0.7492265, -0.2482998, 0.251698, 0.750771, 2.74560E-6, 4.07323E-2, 
      -0.7489055, -0.2475950, 0.252402, 0.751089, 5.49121E-6, 5.76037E-2, 
      -0.7484509, -0.2465978, 0.253396, 0.751539, 1.09824E-5, 8.14629E-2, 
      -0.7478068, -0.2451864, 0.254801, 0.752174, 2.19648E-5, 0.115202, -0.7468933, 
      -0.2431882, 0.256787, 0.753069, 4.39296E-5, 0.162912, -0.7455965, -0.2403575, 
      0.259594, 0.754329, 8.78593E-5, 0.230367, -0.7437523, -0.2363443, 0.263560, 
      0.756099, 1.75718E-4, 0.325717, -0.7411238, -0.2306471, 0.269151, 0.758589, 
      3.51437E-4, 0.460432, -0.7373473, -0.2225595, 0.277014, 0.762103, 7.02874E-4, 
      0.650580, -0.7319174, -0.2110147, 0.288131, 0.766983, 1.40575E-3, 0.918450, 
      -0.7240245, -0.1944645, 0.303826, 0.773772, 2.81149E-3, 1.294163, -0.7123918, 
      -0.1705464, 0.326018, 0.783146, 4.21724E-3, 1.579253, -0.7030419, -0.1516932, 
      0.343154, 0.790253, 4.92012E-3, 1.702698, -0.6988340, -0.1433086, 0.350682, 
      0.793341, 5.27156E-3, 1.760867, -0.6968153, -0.1393073, 0.354255, 0.794798, 
      5.44728E-3, 1.789165, -0.6958245, -0.1373485, 0.355999, 0.795508, 5.53513E-3, 
      1.803129, -0.6953335, -0.1363790, 0.356861, 0.795858, 5.57906E-3, 1.810066, 
      -0.6950891, -0.1358967, 0.357290, 0.796033, 5.60103E-3, 1.813523, -0.6949671, 
      -0.1356561, 0.357504, 0.796119, 5.61201E-3, 1.815249, -0.6949062, -0.1355359, 
      0.357611, 0.796163, 5.61750E-3, 1.816111, -0.6948758, -0.1354759, 0.357664, 
      0.796184, 5.62025E-3, 1.816542, -0.6948606, -0.1354459, 0.357691, 0.796195, 
      5.62162E-3, 1.816758, -0.6948530, -0.1354309, 0.357704, 0.796201, 5.62231E-3, 
      1.816866, -0.6948492, -0.1354234, 0.357711, 0.796203, 5.62265E-3, 1.816919, 
      -0.6948473, -0.1354196, 0.357714, 0.796205, 5.62282E-3, 1.816946, -0.6948463, 
      -0.1354177, 0.357716, 0.796205, 5.62291E-3, 1.816960, -0.6948458, -0.1354168, 
      0.357717, 0.796206, 5.62295E-3, 1.816967, -0.6948456, -0.1354163, 0.357717, 
      0.796206, 5.62297E-3, 1.816970, -0.6948455, -0.1354161, 0.357717, 0.796206, 
      5.62298E-3, 1.816972, -0.6948454, -0.1354160, 0.357717, 0.796206, 5.62299E-3, 
      1.816972, -0.6948454, -0.1354159, 0.357717, 0.796206, 5.62299E-3, 1.389091, 
     -0.7730128, -0.3296001, 0.141945, 0.683462, 1.96520E-2, 3.031403, -0.6584261, 
     -0.1438449, 0.249650, 0.704010, 3.36810E-2, 3.992577, -0.5952508, -0.0201336, 
      0.330225, 0.723319, 3.623176E-2, 4.104700, -0.5878491, -7.64257E-3, 0.334665, 
      0.722183, 6.428978E-2, 5.109331, -0.4935578, 0.136658, 0.388008, 0.685795, 
      6.684051E-2, 5.443877, -0.4637176, 0.169525, 0.408346, 0.686020, 0.080869, 
      9.192625, -0.1945170, 0.376827, 0.601918, 0.830342, 0.087884, 13.26475, 
     -0.024133, 0.503454, 0.716845, 0.893523, 0.094898, 20.72019, 0.158521, 
      0.634920,0.810799,0.942288, 0.097449, 25.56530, 0.243899, 0.717475, 
      0.876807,0.967969,0.098724, 26.92882, 0.268130, 0.741405, 0.896438, 0.978283, 
      0.099362,27.69197,0.280996, 0.753823, 0.906756, 0.983387, 0.099681, 28.08663, 
      0.287506,0.760056,0.911978, 0.985863, 0.099840, 28.28787, 0.290783, 0.763185, 
      0.914612,0.987079,0.099920, 28.37901, 0.292300, 0.764642, 0.915840, 0.987638, 
      0.099960,28.40346,0.292797, 0.765143, 0.916259, 0.987826, 0.100000, 34.40220, 
      0.305826,0.793149,0.929877, 0.986042, 0.100039, 28.36011, 0.292929, 0.765222,
      0.916287,0.987830,0.100079, 28.29237, 0.292567, 0.764802, 0.915895, 0.987646, 
      0.100159,28.11547,0.291319, 0.763508, 0.914726, 0.987097, 0.100318, 27.74551, 
      0.288576,0.760710,0.912212, 0.985901, 0.100637, 27.02355, 0.283149, 0.755169,
      0.907255,0.983479,0.101275, 25.64506, 0.272472, 0.744241, 0.897557, 0.978525, 
      0.102550,23.19311,0.252783, 0.723771, 0.879604, 0.968727, 0.105101, 17.67882, 
      0.191279,0.670074,0.832921, 0.945668, 0.112115, 8.391549, 0.0124316,0.551601, 
      0.745790,0.904137,0.119130, 4.525436,-0.149199, 0.466066, 0.659391, 0.854440, 
      0.126144,2.670650,-0.291271, 0.399669, 0.593693, 0.803774, 0.133159,1.665226, 
     -0.394152,0.312729,0.545130, 0.759130, 0.135710, 1.446525,-0.415016, 0.288925, 
      0.539501,0.762389,0.149739, 0.809176,-0.441209, 0.255551, 0.559359, 0.797353,
      0.163768,0.472591,-0.447064, 0.238353, 0.569192, 0.813019, 0.166318,0.428708, 
     -0.447854,0.236636,0.570986, 0.814478, 0.180347, 0.233970,-0.479907, 0.179081, 
      0.543565,0.803918,0.194377, 0, 0, 0, 0, 0 };
      correctVals = {  54.11599, 0.44944, 9.3508E-2, 2.0747E-2 };
      REQUIRE( ranges::equal( correctSCR, scr, equal ) );
      REQUIRE( ranges::equal( correctVals, out, equal ) );
    } // WHEN

    WHEN( "High tolerance and high energy " ){
      double enow = 0.2;
      tol = 1e-1;
      jbeta = -7;
      for ( auto& val : scr ){ val = 0.0; }

      std::vector<double> xsi(100,0.0);
      auto out = prepareEpMu(enow,j,tev,tol,lat,iinc,lasym,alphas,betas,sab,az,boundXsVec,teff,nbin,jbeta,scr,lastVals,y);

      correctSCR = {0, 0, 0, 0, 0, 0, 8.05839E-7, 0.0131679, -0.7495806, -0.249078, 
      0.250920, 0.750418, 1.61167E-6, 0.0186222, -0.7494065, -0.2486961, 0.251302, 
      0.750590, 3.22335E-6, 0.0263356, -0.7491600, -0.2481554, 0.251840, 0.750834, 
      6.44671E-6, 0.0372439, -0.7488106, -0.2473902, 0.252602, 0.751178, 1.28934E-5, 
      0.0526699, -0.748315, -0.2463068, 0.253678, 0.751664, 2.57868E-5, 0.074483,
     -0.747611, -0.244772, 0.255197, 0.752347, 5.15737E-5, 0.105328, -0.746609,
     -0.242596, 0.257343, .753307, 1.03147E-4, 0.148935, -0.7451816, -0.239509, 
      0.260370, 0.754651, 2.06294E-4, 0.210565, -0.7431378, -0.235122, 0.264640, 
      0.756527, 4.12589E-4, 0.297612, -0.7401984, -0.228870, 0.270629, 0.759158, 
      8.25179E-4, 0.420400, -0.7359181, -0.219951, 0.279049, 0.762795, 1.65035E-3, 
      0.593182, -0.7296389, -0.207080, 0.290930, 0.767802, 3.30071E-3, 0.835018, 
     -0.720274, -0.188411, 0.307594, 0.774557, 6.60143E-3, 1.169834, -0.705983, 
     -0.160986, 0.330952, 0.783594, 1.32028E-2, 1.625559, -0.6829625, -0.118649, 
      0.365656, 0.796264, 2.64057E-2, 2.209573, -0.6458400, -5.568446E-2, 0.412059, 
      0.811785, 5.28115E-2, 2.937096, -0.5740339, 5.864191E-2, 0.497444, 0.839564, 
      0.0792172, 3.408197, -0.500947, 0.164217, 0.572751, 0.866473, 0.0924201, 
      3.591839, -0.463113, 0.215092, 0.607220, 0.878551, 9.90215E-2, 3.647072, 
     -0.446950, 0.235032, 0.625850, 0.886813, 0.102322, 3.683744, -0.437510, 
      0.247073, 0.634352, 0.889953, 0.103972, 3.701664, -0.432761, 0.253085, 
      0.638592, 0.891521, 0.104797, 3.710526, -0.430379, 0.256090, 0.640709, 
      0.892305, 0.105210, 3.714934, -0.429186, 0.257592, 0.641767, 0.892697, 
      0.105416, 3.717132, -0.428589, 0.258342, 0.642296, 0.892893, 0.105519, 
      3.718229, -0.428290, 0.258718, 0.642560, 0.892991, 0.105571, 3.718778, 
     -0.428141, 0.258906, 0.642692, 0.893040, 0.105597, 3.719052, -0.428066, 
      0.258999, 0.642758, 0.893065, 0.105610, 3.719189, -0.428029, 0.259046, 
      0.642791, 0.893077, 0.105616, 3.719258, -0.428010, 0.259070, 0.642808, 
      0.893083, 0.105619, 3.719292, -0.428001, 0.259082, 0.642816, 0.893086, 
      0.105621, 3.719309, -0.427996, 0.259087, 0.642820, 0.893088, 0.105622, 
      3.719317, -0.427994, 0.259090, 0.642822, 0.893089, 0.105622, 3.719322, 
     -0.427993, 0.259092, 0.642823, 0.893089, 0.105622, 3.719324, -0.427992, 
      0.259093, 0.642824, 0.893089, 0.105622, 3.719325, -0.427992, 0.259093, 
      0.642824, 0.893089, 0.105622, 3.719325, -0.427992, 0.259093, 0.642824, 
      0.893089, 0.105622, 3.719326, -0.427992, 0.259093, 0.642824, 0.893089, 
      0.105622, 3.719326, -0.427992, 0.259093, 0.642824, 0.893089, 0.105623, 
      2.973236, -0.517234, 0.104622, 0.473816, 0.766176, 0.133681, 3.682241, 
     -0.378941, 0.315937, 0.637128, 0.842343, 0.136231, 3.688598, -0.372817, 
      0.323154, 0.640697, 0.843247, 0.164289, 3.945966, -0.277831, 0.437559, 
      0.679771, 0.834652, 0.166840, 4.109539, -0.252372, 0.467508, 0.690227, 
      0.836427, 0.180869, 6.221135, -4.941113E-3, 0.659307, 0.78934, 0.910889, 
      0.187884, 8.609271, 0.182882, 0.732273, 0.850863, 0.9445554, 0.194898, 
      13.43754, 0.382317, 0.814861, 0.908577, 0.9705781, 0.19744927, 15.97323, 
      0.442959, 0.847824, 0.934566, 0.983121, 0.19872464, 16.72848, 0.461579, 
      0.859629, 0.944732, 0.988585, 0.199362, 17.16709, 0.4716975, 0.865960, 
      0.950088, 0.991220, 0.199681, 17.39746, 0.4768549, 0.869176, 0.952790, 
      0.992492, 0.199840, 17.51400, 0.479428, 0.870781, 0.954136,  0.993111, 
      0.199920, 17.55609, 0.480412, 0.871423, 0.954675, 0.9933547, 0.199960, 
      17.56363, 0.480656, 0.871612, 0.954835, 0.9934268, 0.19998008, 18.38137, 
      0.469944, 0.860670, 0.952171, 0.993069, 0.19999004, 22.07671, 0.525891, 
      0.887747, 0.953547, 0.991187, 0.199992, 22.12073, 0.5264864, 0.888049, 
      0.953780, 0.991354, 0.199993, 22.13113, 0.5266269, 0.888120, 0.953835, 
      0.991393, 0.199993, 22.13607, 0.526693, 0.888154, 0.953861, 0.991411, 
      0.199993, 22.13732, 0.526710, 0.888162, 0.953868, 0.9914164, 0.199993, 
      22.13763, 0.526714, 0.888165, 0.953869, 0.9914176, 0.19999358, 28.59911, 
      0.588375, 0.870643, 0.937584, 0.986687, 0.19999359, 28.59926, 0.588376, 
      0.870643, 0.937585, 0.986688, 0.199993, 28.59988, 0.5883818, 0.870646, 
      0.937588, 0.986690, 0.199993, 28.60217, 0.5884012, 0.870656, 0.937600, 
      0.986700, 0.199995, 28.61963, 0.588549, 0.870731, 0.937688, 0.986775, 
      0.199997, 28.64467, 0.588762, 0.870841, 0.937817, 0.9868832, 0.199998, 
      28.65102, 0.588817, 0.870870, 0.937852, 0.986911, 0.19999907, 28.65188, 
      0.588825, 0.870874, 0.937857, 0.986915, 0.19999911, 28.65197, 0.588826, 
      0.870874, 0.937857, 0.986916, 0.199999, 44.18280, 0.5683730, 0.797216, 
      0.896646, 0.974788, 0.199999, 44.18282, 0.5683731, 0.797216, 0.896646, 
      0.974788, 0.199999, 44.18286, 0.568373, 0.797217, 0.896646, 0.974788, 
      0.199999, 44.18299, 0.568373, 0.797217, 0.896647, 0.974789, 0.199999, 
      44.18319, 0.568374, 0.797218, 0.896649, 0.974791, 0.20000000, 44.18330, 
      0.568376, 0.797221, 0.896652, 0.974794, 0.20000063, 44.18216, 0.568375, 
      0.797218, 0.896649, 0.974791, 0.200000, 44.18168, 0.5683741, 0.797217, 
      0.896647, 0.974789, 0.200000, 44.18141, 0.5683734, 0.797216, 0.896646, 
      0.974788, 0.200000, 44.18137, 0.5683734, 0.797216, 0.896646, 0.974788, 
      0.200000, 28.65102, 0.588826, 0.870874, 0.937857, 0.986916, 0.200000, 
      28.65095, 0.588826, 0.870874, 0.937857, 0.986916, 0.20000094, 28.65085, 
      0.588825, 0.870874, 0.937857, 0.986915, 0.20000125, 28.64964, 0.588818, 
      0.870870, 0.937852, 0.986911, 0.200002, 28.64192, 0.5887636, 0.870841, 
      0.937817, 0.986882, 0.200004, 28.61415, 0.5885521, 0.870731, 0.937688, 
      0.986774, 0.200006, 28.59518, 0.588403, 0.870655, 0.937598, 0.986699, 
      0.200006, 28.59256, 0.588383, 0.870645, 0.937586, 0.986689, 0.200006, 
      28.59223, 0.588380, 0.870644, 0.937584, 0.986687, 0.20000643, 22.13222, 
      0.526720, 0.888166, 0.953870, 0.991417, 0.20000644, 22.13206, 0.526718, 
      0.888165, 0.953869, 0.991417, 0.200006, 22.13140, 0.5267103, 0.888161, 
      0.953866, 0.991414, 0.200006, 22.13024, 0.5266955, 0.888154, 0.953860, 
      0.991410, 0.200006, 22.12504, 0.526629, 0.888120, 0.953834, 0.991392, 
      0.200007, 22.11428, 0.526491, 0.888050, 0.953780, 0.991353, 0.200009, 
      22.06816, 0.525898, 0.887749, 0.953547, 0.991187, 0.20001993, 18.36728, 
      0.469961, 0.860677, 0.952173, 0.993069, 0.20003986, 17.53633, 0.480707, 
      0.871635, 0.954843, 0.993428, 0.200079, 17.50155, 0.4805144, 0.871469, 
      0.954691, 0.993357, 0.200159, 17.40542, 0.4796312, 0.870873, 0.954169, 
      0.993116, 0.200318, 17.18246, 0.477261, 0.869362, 0.952859, 0.992503, 
      0.200637, 16.74548, 0.472514, 0.866339, 0.950233, 0.991247, 0.201275, 
      15.91687, 0.463228, 0.860423, 0.945053, 0.988654, 0.20255073, 14.46048, 
      0.446313, 0.849542, 0.935322, 0.983328, 0.20510146, 11.03339, 0.390139, 
      0.819338, 0.910813, 0.971236, 0.212115, 5.380518, 0.2019028, 0.746750, 
      0.859307, 0.947590, 0.219130, 2.953821, 0.0315764, 0.684927, 0.806993, 
      0.918380, 0.226144, 1.798367, -0.089644, 0.610489, 0.765750, 0.887029, 
      0.233159, 1.148179, -0.181371, 0.521886, 0.732424, 0.859167, 0.235710, 
      1.001930, -0.200881, 0.498213, 0.726343, 0.859025, 0.24973923, 0.562796, 
     -0.225543, 0.462190, 0.725229, 0.873616, 0.26376824, 0.3273604, -0.232957, 
      0.447126, 0.721330, 0.879156, 0.266318, 0.2601356, -0.2925191, 0.368554, 
      0.674735, 0.883254, 0.280347, 0.149717, -0.3029577, 0.350598, 0.658237, 
      0.871978, 0.294377, 0.0, 0.0, 0.0, 0.0, 0.0 };
      correctVals = { 45.34451, 0.48183, 0.17797, 0.10008 };
      REQUIRE( ranges::equal( correctVals, out, equal ) );
      REQUIRE( ranges::equal( correctSCR, scr, equal ) );

    } // WHEN
  } // GIVEN
} // TEST CASE


